Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Inconsistencies and improvements to SST model #2329

Open
wants to merge 52 commits into
base: develop
Choose a base branch
from

Conversation

rois1995
Copy link
Contributor

@rois1995 rois1995 commented Jul 25, 2024

Proposed Changes

Hi everyone,

I have found some inconsistencies with respect to the literature on the implementation of the Menter's SST model. I would like to use this branch as test bench for any corrections/improvements made to the SST model.

Implementation errors found:

  • Use of production limiter constant in the clipping of the cross-diffusion term for the computation of the F1 blending function in CTurbSSTVariable.cpp.
  • SST naming conventions (V in V2003m should stay for Vorticity production term.). Not yet assessed.
  • Wrong cross-diffusion term included into the residual of w in turb_sources.hpp. It should not be only the positive value as considered in the computation of the F1 blending function. As stated in https://doi.org/10.2514/3.12149. "CDkw is the
    positive portion of the cross-diffusion term in Eq (A2)" pag. 1604. Moreover, a clipping has been introduced for large negative values of this term, as suggested in Peng, Shia-Hui, Peter Eliasson, and Lars Davidson. "Examination of the shear stress transport assumption with a low-Reynolds number k-omega model for aerodynamic flows." Eq 17.
  • Boundary conditions errors at inlets, following the Issue Turbulent Kinetic Energy(TKE) on energy equation in SST model. #1851. A fix has been proposed following @pcarruscag suggestions for the supersonic inlet BC.
  • Boundary conditions definitions are different than the ones proposed in TMR. Maybe give the user the choice to use the BCs implemented in SU2 or the ones from TMR.

Changes to SST model proposed:

  • Inclusion of non modified versions of SST. In this case I think that more work is needed to be sure that the correct terms are taken into account everywhere in the code.
  • Give the user the possibility of changing the production limiter for TKE (not essential).
  • Change values of default turbulent to laminar viscosity to 1e-2 (NASA TMR reports that it should be in the range of 1e-2 to 1e-5, but in SU2 this is set to 10 as default).

I've seen the proposed changes to the lower limits of k and w in #2323 and I tried implementing it in my branch. However, if the implementation proposed in the respective PR is preferred then I will change mine.

Related Work

#2323 #1851

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

  • I am submitting my contribution to the develop branch.
  • My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
  • My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
  • I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
  • I have added a test case that demonstrates my contribution, if necessary.
  • I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

@rois1995 rois1995 changed the title Inconsistencies and improvements to SST model [WIP] Inconsistencies and improvements to SST model Jul 25, 2024
Comment on lines 859 to 860
nPrandtl_Lam, /*!< \brief Number of species
addDoubleOption("FREESTREAM_TURB2LAMVISCRATIO", TurbIntensityAndViscRatioFreeStream[1], 10.0); Prandtl number. */
Copy link
Contributor

@bigfooted bigfooted Jul 25, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Was this intentionally commented, or...? If it is used, I guess this should go to line 872

Suggested change
nPrandtl_Lam, /*!< \brief Number of species
addDoubleOption("FREESTREAM_TURB2LAMVISCRATIO", TurbIntensityAndViscRatioFreeStream[1], 10.0); Prandtl number. */
nPrandtl_Lam, /*!< \brief Number of species laminar Prandtl number. */
addDoubleOption("FREESTREAM_TURB2LAMVISCRATIO", TurbIntensityAndViscRatioFreeStream[1], 10.0); /*!<\brief Freestream mu_turb to mu_lam viscosity ratio */

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is just a wrong copy-paste from me. It should not be there in the first place.

@bigfooted
Copy link
Contributor

Great contribution, Thanks @rois1995 !

@pcarruscag
Copy link
Member

If you are looking into robustness aspects too you should get in touch with @emaberman and @YairMO, seems like they have some good ideas and between the free time of 3 people a lot more can get done :)

@YairMO
Copy link

YairMO commented Jul 25, 2024

Hi,

Regarding the cross-diffusion term (CD) that appears in Omega source term (residual). The SST model (1994/2003) is a high-Reynolds-number model. Namely, It can not predict correctly the sub-layer region (especially the correct profile of the TKE). Therefore, only a positive contribution is required. Moreover, since the SST model was design as a k-w and k-epsilon blending, the CD term "belongs" only to the k-epsilon "branch", that is why the CD term include the factor "1-F1". However, it may happen, that the factor "1-F1" is not a 100% safe guarantee. It may happen that "1-F1" is not zero in region where the CD term is negative (this happen due numerical errors). To avoid such a situation, it is a good idea to clip the CD term with zero. Otherwise, severe numerical robustness issues may rise.
Yes, it is different from Menter publications, but I think that clipping the CD term with zero is completely inline with Menter original idea (that is why, I think, he was including the factor "1-F1". But again the 1-F1 factor is not 100% percent "safe").

- Given option for cross diffusion limiting in W residual
@YairMO
Copy link

YairMO commented Jul 26, 2024

Hi,

The use of an Omega production limiter (about the cross-diffusion term) is correct for low-Reynolds-number (LRN) models (the approach described by Peng et al. is very naive; there are other more rigorous treatments). For high-Reynolds-number (HRN) models, the clipping should be zero, keeping the cross-diffusion term positive; thus, the current implementation is correct.

Indeed, it is not exactly as it appears in Menter's original publication. The factor (1-F1) aimed to promise that the cross-diffusion term will be activated only outside the boundary layer, where it is positive (the cross-diffusion term switches its sign deep inside the boundary layer). This was also recognized by Peng et al. (first paragraph above Eq. 17). However, it may happen that the factor (1-F1)=1 where the cross-diffusion term is negative. Usually, it may happen at the wake, very near the airfoil trailing edge, where the upper and lower boundary layers merge. It is due to the imperfection of the F1 function.

To summarize, the current implementation is correct, and it is perfect for HRN models.

@YairMO
Copy link

YairMO commented Jul 26, 2024

For the sake of clarity, "current implementation" refers to the current treatment of the production code

@emaberman
Copy link
Contributor

What YairMO is saying, is that allowing negative cross diffusion values is incorrect for high Reynolds models and should not be an option, this is a fix used for low Reynolds models only

@rois1995
Copy link
Contributor Author

rois1995 commented Jul 26, 2024

Hi @YairMO, Hi @emaberman ,

thank you very much for your comments. I haven't found any suggestion in literature to clip to only positive values the cross-diffusion term in the w-equation. I understand that it might be more robust, but it is not the standard implementation of the SST model, which is the first thing that we need to achieve. Only then we can build on top of that to improve the robustness of SU2.

Nevertheless, I tried the SWBLI test case and I compared the results across 6 different combinations:

1- develop branch, no changes
2- develop branch, changes to Supersonic_inlet profile as suggested in #1851
3- my branch, with original CDkw implementation (should give exactly the same result as develop+modified BC)
3- my branch, with original CDkw implementation and using boundary conditions from TMR
4- my branch, with original CDkw implementation and using your suggestions for lower limits for k and w.
5- my branch, allowing negative values of CDkw
6- my branch, allowing negative values of CDkw and using boundary conditions from TMR
7- my branch, allowing negative values of CDkw and using your suggestions for lower limits for k and w.
8- my branch, allowing negative values of CDkw, using boundary conditions from TMR and using your suggestions for lower limits for k and w.

When my branch is used, then the changes to the supersonic inlet BC are already in place.

I haven't achieved convergence with 1, 2 and 3. More precisely, 1 diverged right away (after 30 iterations), while 2 and 3 gave "FGMRES - Orthogonalization Failed" after 900ish iterations.

Here you can see the residuals for the different combinations.

OrigCDkw

NegCDkw

Unfortunately I will be busy with the AIAA Conference next week, thus I don't know how much I will be able to work on this. The next test case will be the 2D airfoil near-wake from TMR.

@YairMO
Copy link

YairMO commented Jul 26, 2024

Hi rois1995,

First of all, enjoy your time in Las Vegas. Any paper that you are presenting?

As for our discussion about the cross-diffusion term, I've emailed the "source" (Menter). I believe he will make it clear.
It may be that he will be able to answer only in a while ...

@rois1995
Copy link
Contributor Author

rois1995 commented Sep 28, 2024

Sure, thank you for checking!

Regarding the SWBLI case, if I understood correctly, you have used your code with NK solver from the start and you changed the number of linear solver iterations, right? Can you also try with SST or is your implementation different only for SA?

@YairMO
Copy link

YairMO commented Sep 28, 2024

Yes, your understanding is correct.
As for the SST ... checking

@YairMO
Copy link

YairMO commented Sep 28, 2024

Here are the results obtained using the SST2003 model

SWBLI-Conv-SST2003

@rois1995
Copy link
Contributor Author

It seems that yoru implementation affects also the SST model, is it correct? Nevertheless, huge improvement wrt both this branch and develop branch results (density residual stagnates around -10)!

@YairMO
Copy link

YairMO commented Sep 28, 2024

Yes, it also affects the SST model. Thanks

@YairMO
Copy link

YairMO commented Sep 29, 2024

Hi Rois1995,

I have permission to share the VFE-2 grid; please send me your private email.

@rois1995
Copy link
Contributor Author

That's perfect, thank you! [email protected]

@rois1995
Copy link
Contributor Author

rois1995 commented Oct 4, 2024

I have a few updates. I found out that most of my FGMRES problems in this branch where related to the fact that I was using single precision for the linear system. You can see it in this graphs for the residuals of Rho.

RMSRho_Mesh_3

The first two curves are in mixed-precision and they stop long before reaching the minimum residual required due to divergence of the linear solver.

Moreover, it seems that the 2003 model here implemented (which, simply put, considers the full reynolds stress tensor for the computation of the production of k and always considers k in the stress tensor and in the thermodynamic variables) has faster convergence and increases the recirculating zone (maybe due to the reduction of turbulence kinetic energy) which is a feature that has also been seen in [DOI:10.1017/aer.2020.93].

SFC_Mesh_3

There are some differences between the use of the TMR boundary conditions and the ones used before.

The results improve with respect to the develop in the prediction of the SFC distribution, which has results completely off from the V&V page of SU2

SFC_Mesh_3_Old

RMSRho_Mesh_3_Old

I may keep on checking if the implementation of the v2003 model is correct or not by searching for other test cases (probably coming from DOI:10.1017/aer.2020.93).

@YairMO
Copy link

YairMO commented Oct 4, 2024

How can one switch between single or double precision?

@rois1995
Copy link
Contributor Author

rois1995 commented Oct 4, 2024

when you execute the meson command use the flag -Denable-mixedprec=true

@bigfooted
Copy link
Contributor

bigfooted commented Oct 5, 2024

I am trying to catch up on the developments in this PR but I've lost the overview. Are there still issues that we need to figure out? It looks like convergence is now recovered and results look OK?
Regarding the symmetry plane, as @pcarruscag mentioned it is improving the accuracy, but not necessarily the convergence. @pcarruscag and I tried to get it to converge as well as before, but no guarantees. Do you still see issues with the new symmetry plane?
Regarding the peak residual that you observed in a single node in the mesh: This is a common observation when using Green-Gauss, which gives O(1) accuracy in mixed mesh nodes where quadrilateral and triangular meshes meet. There are some methods that try to fix this, although none of them are holy grails and the best methods are the most expensive ones.

@rois1995
Copy link
Contributor Author

rois1995 commented Oct 7, 2024

The 2003 version of SST seems to be working fine. However, the 1994 version is not. I get Mach 7 values on the upper symmetry boundary even on the coarser mesh. I am trying to understand what is happening.

@rois1995
Copy link
Contributor Author

Maybe a stupid question, but I noticed that if I print the coordinates of the i-th point in the CTurbSSTSolver class and the problem is 2D, then the z-coordinate is not 0 but has random values. Does this affects anything in the code or it is not important since the problem is 2D?

@pcarruscag
Copy link
Member

pcarruscag commented Oct 10, 2024

The z coordinate is not allocated for 3D, so you are either reading out of bounds or reading the x or y value from somewhere else.

@rois1995
Copy link
Contributor Author

That's what I thought, thanks for confirming it.

@rois1995
Copy link
Contributor Author

One more thing: should I change the naming scheme for the model version of SST? As it is now, the naming scheme is V+model version but the NASA TMR site says to use the V only for the vorticity form of the production term. Should we remove the V then, or maybe set it as lower case? This will break all of the configs, no problem for me to go through them all and change it though.

@YairMO
Copy link

YairMO commented Oct 11, 2024

Chris Rumsey made a tremendous effort to standardize popular RANS turbulence models (via the TMR website).
In my opinion, follow his work. Namely, change the present names to exactly what is given on the TMR.

@pcarruscag
Copy link
Member

No changes that break backward compatibility.

@rois1995
Copy link
Contributor Author

That's what I thought.

@@ -865,7 +865,7 @@ class CConfig {
nMu_Temperature_Ref, /*!< \brief Number of species reference temperature for Sutherland model. */
nMu_S, /*!< \brief Number of species reference S for Sutherland model. */
nThermal_Conductivity_Constant,/*!< \brief Number of species constant thermal conductivity. */
nPrandtl_Lam, /*!< \brief Number of species laminar Prandtl number. */
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The original description was correct, for species transport we can have a Prandtl number for each of the species

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah I do not recall changing that but for sure it was wrong in my branch. Corrected, thanks!

su2double Omega_Freestream = 10 * ModVel_FreeStream / config->GetLDomain();
Tke_FreeStream = Omega_Freestream*(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream())/Density_FreeStream;
} else if (config->GetSSTParsedOptions().sust) {
su2double Omega_Freestream = 5*ModVel_FreeStream/config->GetLength_Reynolds();

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable Omega_Freestream is not used.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants